This is the same code as the original analysis, but I’m now modelling the raw amount of buckthorn (instead of the relative abundance)
library(plyr)
library(tidyverse)
library(nlme)
library(lme4)
library(glmmTMB)
library(RColorBrewer)
library(ggrepel)
library(gridExtra)
library(DHARMa)
library(emmeans)
library(dplyr)
gap.area <- read.csv("data/gap_area.csv") %>% select(-X, -Gap)
names(gap.area)[4] <-"Gap Area"
str(gap.area)
## 'data.frame': 48 obs. of 4 variables:
## $ Location: chr "Code Farm" "Code Farm" "Code Farm" "Code Farm" ...
## $ Category: chr "EAB Gap" "EAB Gap" "EAB Gap" "EAB Gap" ...
## $ ID : int 1 2 3 4 1 2 3 4 1 2 ...
## $ Gap Area: num 90.4 139.6 114.1 166.4 59.8 ...
dist.edge <- read.csv("data/dist_edge.csv")
names(dist.edge)[4]<-"Edge Distance"
str(dist.edge)
## 'data.frame': 72 obs. of 4 variables:
## $ Location : chr "Code Farm" "Code Farm" "Code Farm" "Code Farm" ...
## $ Category : chr "EAB Gap" "Other Gap" "Non-Gap" "EAB Gap" ...
## $ ID : int 1 1 1 2 2 2 3 3 3 4 ...
## $ Edge Distance: num 5.59 91.01 85.46 61.89 54.86 ...
frag.size <-read.csv("data/frag_size.csv")
names(frag.size)[2] <- "Fragment Size"
str(frag.size)
## 'data.frame': 6 obs. of 2 variables:
## $ Location : chr "Westminster Ponds" "Meadowlily Woods" "Medway Valley" "Five Points Forest" ...
## $ Fragment Size: int 199 324 230 174 11 37
shade.tol <- read.csv("data/shade_tol.csv") %>% select(-X, -Gap, -tot.BA)
names(shade.tol)[4] <- "Shade Tolerance"
str(shade.tol)
## 'data.frame': 24 obs. of 4 variables:
## $ Location : chr "Code Farm" "Code Farm" "Code Farm" "Code Farm" ...
## $ Category : chr "Non-Gap" "Non-Gap" "Non-Gap" "Non-Gap" ...
## $ ID : int 1 2 3 4 1 2 3 4 1 2 ...
## $ Shade Tolerance: num 4.44 4.52 4.18 4.53 3.33 ...
tree <- read.csv("data/percent_t.csv") %>% select(-X) #wondering if this should be absolute or relative...
names(tree)[4] <- "Percent Tree" #absolute right now (not relative to shrub cover)
str(tree)
## 'data.frame': 72 obs. of 4 variables:
## $ Location : chr "Code Farm" "Code Farm" "Code Farm" "Code Farm" ...
## $ ID : int 1 1 1 2 2 2 3 3 3 4 ...
## $ Category : chr "EAB Gap" "Non-Gap" "Other Gap" "EAB Gap" ...
## $ Percent Tree: int 50 70 60 90 50 40 80 110 120 30 ...
mgmt <- read.csv("data/management.csv")
str(mgmt)
## 'data.frame': 6 obs. of 3 variables:
## $ Location : chr "Westminster Ponds" "Meadowlily Woods" "Medway Valley" "Five Points Forest" ...
## $ Management : chr "Public" "Public" "Public" "Private" ...
## $ Buckthorn_Management: int 1 1 1 1 0 0
EI <- full_join(dist.edge, gap.area, by = c("Location", "Category", "ID"))
EI2 <- full_join(EI, frag.size, by = c("Location"))
EI3 <- full_join(EI2, shade.tol, by = c("Location", "Category", "ID"))
EI4 <- full_join(EI3, tree, by = c("Location", "Category", "ID"))
EI5 <- full_join(EI4, mgmt, by = c("Location"))
EI5$Category <- as.factor(EI5$Category)
names(EI5) <- sub(" ", ".", names(EI5)) #Replace spaces with .
str(EI5)
## 'data.frame': 72 obs. of 10 variables:
## $ Location : chr "Code Farm" "Code Farm" "Code Farm" "Code Farm" ...
## $ Category : Factor w/ 3 levels "EAB Gap","Non-Gap",..: 1 3 2 1 2 3 3 1 2 1 ...
## $ ID : int 1 1 1 2 2 2 3 3 3 4 ...
## $ Edge.Distance : num 5.59 91.01 85.46 61.89 54.86 ...
## $ Gap.Area : num 90.4 59.8 NA 139.6 NA ...
## $ Fragment.Size : int 37 37 37 37 37 37 37 37 37 37 ...
## $ Shade.Tolerance : num NA NA 4.44 NA 4.52 ...
## $ Percent.Tree : int 50 60 70 90 50 40 120 80 110 30 ...
## $ Management : chr "Private" "Private" "Private" "Private" ...
## $ Buckthorn_Management: int 0 0 0 0 0 0 0 0 0 0 ...
Variables:
shade.tol.site <- EI5 %>% #summarize shade tolerance by site
filter(Category == "Non-Gap") %>%
group_by(Location) %>%
summarize(Shade.Tolerance.Site = mean(Shade.Tolerance))
str(shade.tol.site)
## tibble [6 × 2] (S3: tbl_df/tbl/data.frame)
## $ Location : chr [1:6] "Code Farm" "Field Research Station" "Five Points Forest" "Meadowlily Woods" ...
## $ Shade.Tolerance.Site: num [1:6] 4.42 3.33 4.18 4.02 4.82 ...
EI6 <- full_join(EI5, shade.tol.site, by = c("Location")) %>% select(-Shade.Tolerance)
gap.area.site <- EI5 %>% #summarize gap area
filter(Category != "Non-Gap") %>%
group_by(Location) %>%
summarize(Gap.Area.Site = mean(Gap.Area))
EI7 <- full_join(EI6, gap.area.site, by = c("Location")) %>% select(-Gap.Area)
head(EI7)
## Location Category ID Edge.Distance Fragment.Size Percent.Tree Management
## 1 Code Farm EAB Gap 1 5.59 37 50 Private
## 2 Code Farm Other Gap 1 91.01 37 60 Private
## 3 Code Farm Non-Gap 1 85.46 37 70 Private
## 4 Code Farm EAB Gap 2 61.89 37 90 Private
## 5 Code Farm Non-Gap 2 54.86 37 50 Private
## 6 Code Farm Other Gap 2 120.79 37 40 Private
## Buckthorn_Management Shade.Tolerance.Site Gap.Area.Site
## 1 0 4.418994 136.4042
## 2 0 4.418994 136.4042
## 3 0 4.418994 136.4042
## 4 0 4.418994 136.4042
## 5 0 4.418994 136.4042
## 6 0 4.418994 136.4042
#write.csv(EI7, file = "EI.csv")
EI7 %>% select(Edge.Distance, Fragment.Size, Percent.Tree, Shade.Tolerance.Site) -> EI.data
prcomp(EI.data, scale. = TRUE) -> EI.pca
summary(EI.pca)
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.3486 1.0648 0.8749 0.53106
## Proportion of Variance 0.4547 0.2835 0.1914 0.07051
## Cumulative Proportion 0.4547 0.7381 0.9295 1.00000
biplot(EI.pca)
PC1 <- predict(EI.pca)[,1]
PC2 <- predict(EI.pca)[,2]
EI7 %>% select(Edge.Distance, Fragment.Size, Percent.Tree, Shade.Tolerance.Site, Gap.Area.Site) %>%
prcomp(scale. = TRUE) -> EI.pca2
summary(EI.pca2)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.4839 1.1522 0.9228 0.58280 0.52837
## Proportion of Variance 0.4404 0.2655 0.1703 0.06793 0.05584
## Cumulative Proportion 0.4404 0.7059 0.8762 0.94416 1.00000
biplot(EI.pca2)
PC1b <- predict(EI.pca2)[,1]
PC2b <- predict(EI.pca2)[,2]
rescale.U <- rep(EI.pca$sdev, each = length(EI.pca$sdev)) #get lengths
U.scale2 <- EI.pca$rotation * rescale.U #multiply lengths by sqrt SD
round(U.scale2^2,2) #variability in each variable for each PC
## PC1 PC2 PC3 PC4
## Edge.Distance 0.60 0.23 0.07 0.11
## Fragment.Size 0.83 0.00 0.03 0.14
## Percent.Tree 0.03 0.78 0.18 0.02
## Shade.Tolerance.Site 0.37 0.13 0.49 0.01
qqnorm(EI.data[,1]);qqline(EI.data[,1])
qqnorm(EI.data[,2]);qqline(EI.data[,2])
qqnorm(EI.data[,3]);qqline(EI.data[,3])
qqnorm(EI.data[,4]);qqline(EI.data[,4])
#display.brewer.all(colorblindFriendly = TRUE)
mypalette <-brewer.pal(11,"BrBG")[c(1,3,7,5,9,11)]
#mypalette <- c("","","lightseagreen", "lightskyblue2", "dodgerblue", "dodgerblue4")
image(1:9,1,as.matrix(1:9),col=mypalette,ylab="",xaxt="n",yaxt="n",bty="n")
U <- data.frame(EI.pca$rotation)
colnames(U) <- colnames(EI.pca$rotation)
rownames(U) <- rownames(EI.pca$rotation)
U$descriptors <- rownames(U)
F.1 <- data.frame(EI.pca$x) # The book calls this matrix F but I use F.1 because in R, F is shorthand for FALSE
colnames(F.1) <- colnames(EI.pca$x)
rownames(F.1) <- rownames(EI.pca$x)
str(U)
## 'data.frame': 4 obs. of 5 variables:
## $ PC1 : num 0.574 0.674 0.118 0.449
## $ PC2 : num -0.44703 0.00827 0.82678 0.34137
## $ PC3 : num 0.292 0.198 0.486 -0.799
## $ PC4 : num -0.621 0.711 -0.257 -0.207
## $ descriptors: chr "Edge.Distance" "Fragment.Size" "Percent.Tree" "Shade.Tolerance.Site"
levels(U$descriptors) <- c("Distance from Forest Edge", "Forest Fragment Size", "Tree Regeneration", "Successional Stage")
U$descriptors <- as.factor(U$descriptors)
str(U)
## 'data.frame': 4 obs. of 5 variables:
## $ PC1 : num 0.574 0.674 0.118 0.449
## $ PC2 : num -0.44703 0.00827 0.82678 0.34137
## $ PC3 : num 0.292 0.198 0.486 -0.799
## $ PC4 : num -0.621 0.711 -0.257 -0.207
## $ descriptors: Factor w/ 4 levels "Edge.Distance",..: 1 2 3 4
F.1$Location<- EI7$Location
F.1$Category <- EI7$Category
F.1$Location = factor(F.1$Location, levels=c("Field Research Station", "Code Farm", "Westminster Ponds", "Five Points Forest", "Medway Valley", "Meadowlily Woods"))
F.1$Category = factor(F.1$Category, levels=c("EAB Gap", "Other Gap", "Non-Gap"))
#Change Names
F.1$Location <- revalue(F.1$Location, c("Field Research Station"="Private 1"))
F.1$Location <- revalue(F.1$Location, c("Code Farm"="Private 2"))
F.1$Location <- revalue(F.1$Location, c("Westminster Ponds"="Public 1"))
F.1$Location <- revalue(F.1$Location, c("Five Points Forest"="Private 3"))
F.1$Location <- revalue(F.1$Location, c("Medway Valley"="Public 2"))
F.1$Location <- revalue(F.1$Location, c("Meadowlily Woods"="Public 3"))
levels((F.1$Location))
## [1] "Private 1" "Private 2" "Public 1" "Private 3" "Public 2" "Public 3"
Plot the objects
levels(U$descriptors) <- c("Distance from Edge", "Fragment Size", "Tree Regeneration",
"Successional Stage")
biplot1 <- ggplot(F.1, aes(x = PC1, y =PC2)) +
geom_point(aes(shape = Category, fill = Location),col = "black", size = 2, pch=21, alpha = 0.9) +
theme_classic() +
coord_fixed() +
labs(x = 'Principle Component 1', y = "Principle Component 2") +
scale_fill_manual(values = mypalette) +
geom_segment(data = U, aes(xend = PC1*3, yend = PC2*3,x = 0, y = 0), col = "black", alpha = 0.7, arrow = arrow(length = unit(0.35, "cm"))) +
geom_label_repel(data = U, aes(x = PC1*3, y = PC2*3, label = descriptors),
col = "black", nudge_y = -0.35,
segment.colour = NA, size = 3, alpha = 0.7) +
theme(legend.position="bottom", legend.box = "horizontal", plot.margin=grid::unit(c(0,0,0,0), "mm"))
biplot1
#ggsave('figures/fig.biplotBrBu.jpeg',biplot1, units = 'cm', width = 14, height =12)
shrubs.f <- read.csv("data/shrubs_final.csv") %>% select(-"X")
shrubs.EI <- full_join(shrubs.f, EI7, by = c("Location", "Category", "ID"))
## Warning: Column `Category` joining character vector and factor, coercing into
## character vector
shrubs.EI$Category <- as.factor(shrubs.EI$Category)
names(shrubs.EI) <- sub(" ", ".", names(shrubs.EI)) #Replace spaces with .
str(shrubs.EI)
## 'data.frame': 72 obs. of 15 variables:
## $ Location : chr "Code Farm" "Code Farm" "Code Farm" "Code Farm" ...
## $ Category : Factor w/ 3 levels "EAB Gap","Non-Gap",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ ID : int 1 2 3 4 1 2 3 4 1 2 ...
## $ richness : int 2 3 2 3 3 2 3 2 3 1 ...
## $ H : num 0.5 0.965 0.693 1.099 0.736 ...
## $ percentNN : num 0 0 0 0 12.5 ...
## $ percentB : num 0 0 0 0 12.5 ...
## $ Buckthorn : int 0 0 0 0 10 0 0 0 10 0 ...
## $ Edge.Distance : num 5.59 61.89 59.79 14.57 85.46 ...
## $ Fragment.Size : int 37 37 37 37 37 37 37 37 37 37 ...
## $ Percent.Tree : int 50 90 80 30 70 50 110 100 60 40 ...
## $ Management : chr "Private" "Private" "Private" "Private" ...
## $ Buckthorn_Management: int 0 0 0 0 0 0 0 0 0 0 ...
## $ Shade.Tolerance.Site: num 4.42 4.42 4.42 4.42 4.42 ...
## $ Gap.Area.Site : num 136 136 136 136 136 ...
shrubs.EI$PC1 <- PC1
shrubs.EI$PC2 <- PC2
hist(shrubs.EI$PC1)
shrubs.EI %>% summarize(mean = mean(PC1) %>% round(1), sd = sd(PC1) %>% round(1)) #mean already was 0 (centered)
## mean sd
## 1 0 1.3
shrubs.EI$PC1.s <- scale(shrubs.EI$PC1, center = TRUE, scale = TRUE) #standardize
shrubs.EI %>% summarize(mean = mean(PC1.s) %>% round(1), sd = sd(PC1.s))
## mean sd
## 1 0 1
head(shrubs.EI)
## Location Category ID richness H percentNN percentB Buckthorn
## 1 Code Farm EAB Gap 1 2 0.5004024 0.0 0.0 0
## 2 Code Farm EAB Gap 2 3 0.9649629 0.0 0.0 0
## 3 Code Farm EAB Gap 3 2 0.6931472 0.0 0.0 0
## 4 Code Farm EAB Gap 4 3 1.0986123 0.0 0.0 0
## 5 Code Farm Non-Gap 1 3 0.7356219 12.5 12.5 10
## 6 Code Farm Non-Gap 2 2 0.5004024 0.0 0.0 0
## Edge.Distance Fragment.Size Percent.Tree Management Buckthorn_Management
## 1 5.59 37 50 Private 0
## 2 61.89 37 90 Private 0
## 3 59.79 37 80 Private 0
## 4 14.57 37 30 Private 0
## 5 85.46 37 70 Private 0
## 6 54.86 37 50 Private 0
## Shade.Tolerance.Site Gap.Area.Site PC1 PC2 PC1.s
## 1 4.418994 136.4042 -1.3830871 0.12935587 -1.0255878
## 2 4.418994 136.4042 -0.8279143 -0.02767019 -0.6139157
## 3 4.418994 136.4042 -0.8262459 0.24628193 -0.6126785
## 4 4.418994 136.4042 -0.8987823 0.85326388 -0.6664657
## 5 4.418994 136.4042 -1.0833052 -0.10406499 -0.8032933
## 6 4.418994 136.4042 -0.7175934 -0.66407255 -0.5321104
lm.mgmt <- lm(data = shrubs.EI, PC1 ~ Management)
summary(lm.mgmt)
##
## Call:
## lm(formula = PC1 ~ Management, data = shrubs.EI)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.48175 -1.08831 0.07374 0.81699 1.81340
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9127 0.1656 -5.510 5.59e-07 ***
## ManagementPublic 1.8255 0.2343 7.793 4.40e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9939 on 70 degrees of freedom
## Multiple R-squared: 0.4645, Adjusted R-squared: 0.4569
## F-statistic: 60.73 on 1 and 70 DF, p-value: 4.397e-11
anova(lm.mgmt)
## Analysis of Variance Table
##
## Response: PC1
## Df Sum Sq Mean Sq F value Pr(>F)
## Management 1 59.983 59.983 60.728 4.397e-11 ***
## Residuals 70 69.142 0.988
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(lm.mgmt)
## 2.5 % 97.5 %
## (Intercept) -1.243108 -0.5823824
## ManagementPublic 1.358287 2.2926933
ggplot(shrubs.EI, aes(x = Management, y = PC1)) +
geom_boxplot() +
geom_jitter(alpha=0.6, aes(col = Location)) +
theme_classic()
EI.location <- shrubs.EI %>% group_by(Location) %>%
summarize(EI = mean(PC1) %>% round(2))
arrange(EI.location, EI)
## # A tibble: 6 x 2
## Location EI
## <chr> <dbl>
## 1 Field Research Station -2.25
## 2 Code Farm -0.92
## 3 Westminster Ponds -0.01
## 4 Five Points Forest 0.43
## 5 Medway Valley 1.17
## 6 Meadowlily Woods 1.57
EI.cat <- shrubs.EI %>% group_by(Category) %>%
summarize(EI = mean(PC1) %>% round(2))
arrange(EI.cat, EI)
## # A tibble: 3 x 2
## Category EI
## <fct> <dbl>
## 1 EAB Gap -0.03
## 2 Other Gap 0
## 3 Non-Gap 0.03
hist(shrubs.f$Buckthorn) #left skewed
(sum(shrubs.EI$Buckthorn > 0) / nrow(shrubs.EI))*100 # 56% zeros
## [1] 55.55556
#shrubs.EI$Category <- relevel(shrubs.EI$Category, ref = "Non-Gap") #re-level to set non-gap as reference condition - removed
Is it more interesting to consider EAB gaps or non-gaps as the reference category? Decided on EAB gaps
shrubs.B1 <- glmer(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location), family = poisson, data = shrubs.EI)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
E1 <- resid(shrubs.B1, type = "pearson")
N <- nrow(shrubs.EI)
p <- length(fixef(shrubs.B1)) + 1
sum(E1^2) / (N - p)
## [1] 14.59909
Yea, that’s way over 1 (super overdispersed)
shrubs.B2 <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location), family = "nbinom2", data = shrubs.EI)
E3 <- resid(shrubs.B2, type = "pearson")
N <- nrow(shrubs.EI)
p <- length(fixef(shrubs.B2)$cond) + 1 + 1
sum(E3^2) / (N - p)
## [1] 0.5120524
Now likely underdispersed, but better than before
Select ecological integrity (PC1 / PC2) with and without gap area
#PC1 without gap area
shrubsfit.Ba <- glmmTMB(Buckthorn ~ Category * PC1 + (1 | Location), family = "nbinom2", data = shrubs.EI)
#PC1 with gap area
shrubsfit.Bb <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location), family = "nbinom2", data = shrubs.EI)
#PC2 without gap area
shrubsfit.Bc <- glmmTMB(Buckthorn ~ Category * PC2 + (1 | Location), family = "nbinom2", data = shrubs.EI)
#PC2 with gap area
shrubsfit.Bd <- glmmTMB(Buckthorn ~ Category * PC2 + Gap.Area.Site + (1 | Location), family = "nbinom2", data = shrubs.EI)
#gap area in PCA
shrubsfit.Be <- glmmTMB(Buckthorn ~ Category * PC1b + (1 | Location), family = "nbinom2", data = shrubs.EI)
AIC(shrubsfit.Ba, shrubsfit.Bb, shrubsfit.Bc, shrubsfit.Bd, shrubsfit.Be)
## df AIC
## shrubsfit.Ba 8 455.3764
## shrubsfit.Bb 9 453.7808
## shrubsfit.Bc 8 453.2132
## shrubsfit.Bd 9 455.1510
## shrubsfit.Be 8 458.7470
PC1 with gap area (shrubsfit.Bb) is best
shrubsfit.Bb <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location), family = "nbinom2", data = shrubs.EI)
shrubsfit.Bb.no.inter <- glmmTMB(Buckthorn ~ Category + PC1 + Gap.Area.Site + (1 | Location), family = "nbinom2", data = shrubs.EI)
shrubsfit.null <- glmmTMB(Buckthorn ~ 1 + (1 | Location), family = "nbinom2", data = shrubs.EI)
anova(shrubsfit.Bb, shrubsfit.Bb.no.inter)
## Data: shrubs.EI
## Models:
## shrubsfit.Bb.no.inter: Buckthorn ~ Category + PC1 + Gap.Area.Site + (1 | Location), zi=~0, disp=~1
## shrubsfit.Bb: Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df
## shrubsfit.Bb.no.inter 7 455.30 471.23 -220.65 441.30
## shrubsfit.Bb 9 453.78 474.27 -217.89 435.78 5.5155 2
## Pr(>Chisq)
## shrubsfit.Bb.no.inter
## shrubsfit.Bb 0.06344 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(shrubsfit.Bb, shrubsfit.Bb.no.inter, shrubsfit.null)
## df AIC
## shrubsfit.Bb 9 453.7808
## shrubsfit.Bb.no.inter 7 455.2962
## shrubsfit.null 3 463.8500
p-value for interaction is now 0.06 (still marginally significant)
summary(shrubs.B2)
## Family: nbinom2 ( log )
## Formula: Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location)
## Data: shrubs.EI
##
## AIC BIC logLik deviance df.resid
## 453.8 474.3 -217.9 435.8 63
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Location (Intercept) 1.08e-08 0.0001039
## Number of obs: 72, groups: Location, 6
##
## Overdispersion parameter for nbinom2 family (): 0.34
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.07861 2.48643 3.651 0.000261 ***
## CategoryNon-Gap -2.66823 0.58783 -4.539 5.65e-06 ***
## CategoryOther Gap -0.56264 0.50306 -1.118 0.263385
## PC1 -0.73065 0.28456 -2.568 0.010238 *
## Gap.Area.Site -0.04269 0.01725 -2.474 0.013348 *
## CategoryNon-Gap:PC1 -0.49142 0.46269 -1.062 0.288202
## CategoryOther Gap:PC1 0.59645 0.35743 1.669 0.095178 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Effect of gap area now significant
There are two types of fitted values for a GLM - the “link” predictions are estimates of \(\eta\) and the response predictions are estimates of \(\mu\).
shrubs.EI$res <- residuals(shrubsfit.Bb, type = "pearson")
shrubs.EI$fit <- predict(shrubsfit.Bb, type = "response")
ggplot(shrubs.EI,aes(x = fit, y = res)) +
geom_point() + geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
ggplot(shrubs.EI, aes(x = PC1, y = res)) +
geom_point() + geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
simulationOutput1 <- simulateResiduals(fittedModel = shrubsfit.Bb, n = 1000)
plotSimulatedResiduals(simulationOutput = simulationOutput1)
## plotSimulatedResiduals is deprecated, please switch your code to simply using the plot() function
par(mfrow = c(1,3))
names(gap.area) <- sub(" ", ".", names(gap.area)) #Replace spaces with .
gap.test <- lme(Gap.Area~Category,random=~1|Location,data=gap.area)
anova(gap.test)
## numDF denDF F-value p-value
## (Intercept) 1 41 96.89140 <.0001
## Category 1 41 1.28749 0.2631
summary(gap.test)
## Linear mixed-effects model fit by REML
## Data: gap.area
## AIC BIC logLik
## 566.6578 573.9724 -279.3289
##
## Random effects:
## Formula: ~1 | Location
## (Intercept) Residual
## StdDev: 0.004330794 97.93859
##
## Fixed effects: Gap.Area ~ Category
## Value Std.Error DF t-value p-value
## (Intercept) 155.1877 19.99163 41 7.762633 0.0000
## CategoryOther Gap -32.0801 28.27244 41 -1.134678 0.2631
## Correlation:
## (Intr)
## CategoryOther Gap -0.707
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.4588304 -0.6725202 -0.2158717 0.3635209 2.2730709
##
## Number of Observations: 48
## Number of Groups: 6
gap.area %>% group_by(Category) %>% summarize(gap.area = mean(Gap.Area))
## # A tibble: 2 x 2
## Category gap.area
## <chr> <dbl>
## 1 EAB Gap 155.
## 2 Other Gap 123.
intervals(gap.test, which = "fixed")
## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 114.81378 155.1877 195.56161
## CategoryOther Gap -89.17744 -32.0801 25.01724
## attr(,"label")
## [1] "Fixed effects:"
ggplot(gap.area, aes(x = Category, y = Gap.Area)) +
geom_boxplot() +
geom_jitter(alpha=0.6, aes(col = Location)) +
theme_classic()
shrubs.0infl <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location) + Management, ziformula = ~., family = "nbinom2", data = shrubs.EI)
shrubs.hurdle <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location) + Management, ziformula = ~., family = "truncated_nbinom2", data = shrubs.EI)
#summary(shrubs.0infl)
#summary(shrubs.hurdle)
AIC(shrubs.0infl, shrubs.hurdle, shrubs.B2) # Better than neg binomial
## df AIC
## shrubs.0infl 19 427.5488
## shrubs.hurdle 19 427.6976
## shrubs.B2 9 453.7808
Both the zero-inflated and zero-altered model are better than the original negative binomial model
We’re going to use zero-altered because it’s safe to assume that if buckthorn was present, it was recorded
shrubs.hurdle <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location),
ziformula = ~ Category * PC1 + Gap.Area.Site + Management + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
shrubs.hurdle1 <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location),
ziformula = ~ Category * PC1 + Gap.Area.Site + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
AIC(shrubs.hurdle, shrubs.hurdle1)
## df AIC
## shrubs.hurdle 18 428.0959
## shrubs.hurdle1 17 426.9657
Drop management, reduces AIC
426.9657-428.0959
## [1] -1.1302
shrubs.hurdle1a <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location),
ziformula = ~ Category + PC1 + Gap.Area.Site + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
shrubs.hurdle1b <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location),
ziformula = ~ Category * PC1 + Gap.Area.Site + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
AIC(shrubs.hurdle1a, shrubs.hurdle1b)
## df AIC
## shrubs.hurdle1a 15 424.8045
## shrubs.hurdle1b 17 426.9657
424.8045-426.9657
## [1] -2.1612
1a - drop interaction, reduced AIC
shrubs.hurdle2c <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location),
ziformula = ~ Category + PC1 + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
AIC(shrubs.hurdle1a, shrubs.hurdle2c)
## df AIC
## shrubs.hurdle1a 15 424.8045
## shrubs.hurdle2c 14 423.5389
423.5389-424.8045
## [1] -1.2656
2c - drop gap area, reducec AIC
shrubs.hurdle3a <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location),
ziformula = ~ Category + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
AIC(shrubs.hurdle2c, shrubs.hurdle3a)
## df AIC
## shrubs.hurdle2c 14 423.5389
## shrubs.hurdle3a 13 422.0300
422.0300-423.5389
## [1] -1.5089
3b - drop PC1, reduced AIC
shrubs.hurdle4 <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location),
ziformula = ~ 1 + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
AIC(shrubs.hurdle3a, shrubs.hurdle4)
## df AIC
## shrubs.hurdle3a 13 422.0300
## shrubs.hurdle4 11 426.9292
426.9292-422.0300
## [1] 4.8992
Keep Category, dropping increases AIC > 2
shrubs.hurdle5c <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + Management + (1 | Location),
ziformula = ~ Category + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
shrubs.hurdle5d <- glmmTMB(Buckthorn ~ Category * PC1 + Gap.Area.Site + (1 | Location),
ziformula = ~ Category + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
AIC(shrubs.hurdle5c, shrubs.hurdle5d)
## df AIC
## shrubs.hurdle5c 14 421.6317
## shrubs.hurdle5d 13 422.0300
422.0300-421.6317
## [1] 0.3983
Drop management
shrubs.hurdle5 <- glmmTMB(Buckthorn ~ Category * PC1 + (1 | Location),
ziformula = ~ Category + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
AIC(shrubs.hurdle5d, shrubs.hurdle5)
## df AIC
## shrubs.hurdle5d 13 422.0300
## shrubs.hurdle5 12 423.2072
423.2072-422.0300
## [1] 1.1772
Drop gap size
shrubs.hurdle6 <- glmmTMB(Buckthorn ~ Category + PC1 + (1 | Location),
ziformula = ~ Category + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI)
AIC(shrubs.hurdle5, shrubs.hurdle6)
## df AIC
## shrubs.hurdle5 12 423.2072
## shrubs.hurdle6 10 427.9651
427.9651-423.2072
## [1] 4.7579
Keep interaction
sim.m5 <- simulate(shrubs.hurdle5, nsim = 999) %>% t() %>% as.data.frame() %>%
gather(observation, sim_value)
sim.m5$Location <- rep(shrubs.EI$Location, each = 999)
sim.m5$obs.num <- as.numeric(str_extract(sim.m5$observation, "[:digit:]+"))
str(sim.m5)
## 'data.frame': 71928 obs. of 4 variables:
## $ observation: chr "V1" "V1" "V1" "V1" ...
## $ sim_value : num 14 77 0 80 59 36 105 41 0 102 ...
## $ Location : chr "Code Farm" "Code Farm" "Code Farm" "Code Farm" ...
## $ obs.num : num 1 1 1 1 1 1 1 1 1 1 ...
ggplot() +
geom_boxplot(data = sim.m5, aes(x = reorder(observation, obs.num),
y = sim_value + 1,
col = Location), alpha = 0.5) +
scale_y_log10() +
geom_point(data = shrubs.EI, aes(x = 1:72, y = Buckthorn + 1), col = "black") +
theme_classic()
A lot of uncertainty but the real observations are usually in line the middle 50% of the simulations.
A bit more buckthorn than the model predicts at Westminster Ponds but I think that is consistent with previous knowledge.
non.min <- min(subset(shrubs.EI$PC1, shrubs.EI$Category == "Non-Gap"))
non.max<- max(subset(shrubs.EI$PC1, shrubs.EI$Category == "Non-Gap"))
EAB.min <- min(subset(shrubs.EI$PC1, shrubs.EI$Category == "EAB Gap"))
EAB.max <- max(subset(shrubs.EI$PC1, shrubs.EI$Category == "EAB Gap"))
other.min <- min(subset(shrubs.EI$PC1, shrubs.EI$Category == "Other Gap"))
other.max <- max(subset(shrubs.EI$PC1, shrubs.EI$Category == "Other Gap"))
pred <- data.frame(PC1 = c(seq(EAB.min, EAB.max, length.out = 100),
seq(other.min, other.max, length.out = 100),
seq(non.min, non.max, length.out = 100)),
Category = c(rep("EAB Gap", 100), rep("Other Gap", 100), rep("Non-Gap", 100)),
Gap.Area.Site = rep(mean(shrubs.EI$Gap.Area.Site),300), # Added gap area
Management = c(rep("Private", 150), rep("Public", 150)),
Location = rep(NA, 300)) # This NA gives general predictions instead of location specific predictions
pred$fit <- predict(shrubs.hurdle5, newdata = pred, se.fit = TRUE, type = "link")$fit
pred$se <- predict(shrubs.hurdle5, newdata = pred, se.fit = TRUE, type = "link")$se.fit
str(pred)
## 'data.frame': 300 obs. of 7 variables:
## $ PC1 : num -2.28 -2.23 -2.19 -2.14 -2.1 ...
## $ Category : chr "EAB Gap" "EAB Gap" "EAB Gap" "EAB Gap" ...
## $ Gap.Area.Site: num 139 139 139 139 139 ...
## $ Management : chr "Private" "Private" "Private" "Private" ...
## $ Location : logi NA NA NA NA NA NA ...
## $ fit : num 4.29 4.27 4.26 4.24 4.23 ...
## $ se : num 0.404 0.398 0.392 0.386 0.38 ...
pred.upr <- pred$fit + 1.96 * pred$se
pred.lwr <- pred$fit - 1.96 * pred$se
pred$upr <- ifelse(exp(pred.upr) < 100, exp(pred.upr), 100)
pred$lwr <- ifelse(exp(pred.lwr) > 0, exp(pred.lwr), 0)
pred$mean <- exp(pred$fit)
mypalette2 = c(brewer.pal(12, "Paired")[c(12,2,4)])
str(pred)
## 'data.frame': 300 obs. of 10 variables:
## $ PC1 : num -2.28 -2.23 -2.19 -2.14 -2.1 ...
## $ Category : chr "EAB Gap" "EAB Gap" "EAB Gap" "EAB Gap" ...
## $ Gap.Area.Site: num 139 139 139 139 139 ...
## $ Management : chr "Private" "Private" "Private" "Private" ...
## $ Location : logi NA NA NA NA NA NA ...
## $ fit : num 4.29 4.27 4.26 4.24 4.23 ...
## $ se : num 0.404 0.398 0.392 0.386 0.38 ...
## $ upr : num 100 100 100 100 100 100 100 100 100 100 ...
## $ lwr : num 33 32.9 32.8 32.7 32.6 ...
## $ mean : num 72.9 71.8 70.8 69.7 68.6 ...
shrubs.EI$Category = factor(shrubs.EI$Category, levels=c("EAB Gap", "Other Gap", "Non-Gap"))
pred$Category = factor(pred$Category, levels=c("EAB Gap", "Other Gap", "Non-Gap"))
fig.EI <- ggplot(pred, aes(x = PC1, y = mean)) +
geom_ribbon(aes(x = PC1, ymin = lwr, ymax = upr, fill = Category), alpha = 0.15) +
geom_point(data = shrubs.EI, aes(x = PC1, y = Buckthorn, col = Category), alpha = 0.6) +
geom_line(aes(col = Category), size = 1, alpha = 0.8) +
labs(x = "Ecological Integrity (index)",
y = "European buckthorn (% cover)",
col = "Category") +
scale_y_continuous(limits = c(0, 100)) +
scale_colour_manual(values = mypalette2) +
scale_fill_manual(values = mypalette2) +
theme_classic() +
theme(legend.position = "none")
fig.EI
#ggsave('figures/fig.EI.jpeg',fig.EI , units = 'cm', width = 14, height = 10)
fig.EI2 <- ggplot(pred, aes(x = PC1, y = mean)) +
geom_ribbon(aes(x = PC1, ymin = lwr, ymax = upr, fill = Category), alpha = 0.15) +
geom_point(data = shrubs.EI, aes(x = PC1, y = Buckthorn, col = Category), alpha = 0.6) +
geom_line(aes(col = Category), size = 1, alpha = 0.8) +
labs(x = "Ecological Integrity (index)",
y = "European buckthorn (% cover)",
col = "Category") +
scale_y_continuous(limits = c(0, 100)) +
scale_colour_manual(values = mypalette2, name = "Gap Category", labels = c("Emerald Ash Borer Gap", "Other Canopy Gap", "Non-Gap")) +
scale_fill_manual(values = mypalette2, name = "Gap Category", labels = c("Emerald Ash Borer Gap", "Other Canopy Gap", "Non-Gap")) +
theme_classic() +
theme(legend.position = "bottom")
fig.EI2
#ggsave('figures/fig.EI2.jpeg',fig.EI2 , units = 'cm', width = 15, height = 12)
Check those high % buckthorn points - where are they coming from?
ggplot(pred, aes(x = PC1, y = mean)) +
geom_ribbon(aes(x = PC1, ymin = lwr, ymax = upr, fill = Category), alpha = 0.15) +
geom_point(data = shrubs.EI, aes(x = PC1, y = Buckthorn, col = Category, shape = Location), alpha = 0.6) +
geom_line(aes(col = Category), size = 1, alpha = 0.8) +
labs(x = "Ecological Integrity",
y = "Buckthorn (%)",
col = "Category") +
scale_colour_manual(values = mypalette2) +
scale_fill_manual(values = mypalette2) +
theme_classic()
# theme(legend.position = "none")
shrubs.EI.w <- shrubs.EI %>% filter(Location == "Westminster Ponds")
shrubs.EI.nw <- shrubs.EI %>% filter(Location != "Westminster Ponds")
ggplot(pred, aes(x = PC1, y = mean)) +
geom_ribbon(aes(x = PC1, ymin = lwr, ymax = upr, fill = Category), alpha = 0.15) +
geom_point(data = shrubs.EI.nw, aes(x = PC1, y = Buckthorn, col = Category), alpha = 0.6) +
geom_point(data = shrubs.EI.w, aes(x = PC1, y = Buckthorn, col = Category), alpha = 0.6, shape = 17) +
geom_line(aes(col = Category), size = 1, alpha = 0.8) +
labs(x = "Ecological Integrity",
y = "Buckthorn (%)",
col = "Category") +
scale_colour_manual(values = mypalette2) +
scale_fill_manual(values = mypalette2) +
theme_classic()
# theme(legend.position = "none")
The 90% are all coming from Westminster Ponds
summary(shrubs.hurdle5)
## Family: truncated_nbinom2 ( log )
## Formula: Buckthorn ~ Category * PC1 + (1 | Location)
## Zero inflation: ~Category + (1 | Location)
## Data: shrubs.EI
##
## AIC BIC logLik deviance df.resid
## 423.2 450.5 -199.6 399.2 60
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.103 0.3209
## Number of obs: 72, groups: Location, 6
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## Location (Intercept) 1.218 1.104
## Number of obs: 72, groups: Location, 6
##
## Overdispersion parameter for truncated_nbinom2 family (): 2.93
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5180 0.2103 16.728 < 2e-16 ***
## CategoryNon-Gap -1.4585 0.3674 -3.970 7.19e-05 ***
## CategoryOther Gap -0.6244 0.2338 -2.671 0.00756 **
## PC1 -0.3386 0.1579 -2.145 0.03198 *
## CategoryNon-Gap:PC1 0.1869 0.3073 0.608 0.54301
## CategoryOther Gap:PC1 0.5806 0.1908 3.044 0.00234 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.883e-01 6.716e-01 -1.323 0.1859
## CategoryNon-Gap 1.751e+00 7.193e-01 2.435 0.0149 *
## CategoryOther Gap 3.277e-06 6.794e-01 0.000 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(fixef(shrubs.hurdle5)$cond) # Intercept expected cover when present, others are multiplicative changes in cover conditional on presence
## (Intercept) CategoryNon-Gap CategoryOther Gap
## 33.7154190 0.2325890 0.5355936
## PC1 CategoryNon-Gap:PC1 CategoryOther Gap:PC1
## 0.7128030 1.2055000 1.7871821
exp(fixef(shrubs.hurdle5)$zi) # Intercept is odds of zero, others are odds ratios
## (Intercept) CategoryNon-Gap CategoryOther Gap
## 0.4113465 5.7618230 1.0000033
exp(confint(shrubs.hurdle5))
## 2.5 % 97.5 % Estimate
## cond.(Intercept) 22.3264815 50.9139553 33.7154190
## cond.CategoryNon-Gap 0.1132062 0.4778683 0.2325890
## cond.CategoryOther Gap 0.3387368 0.8468538 0.5355936
## cond.PC1 0.5231179 0.9712689 0.7128030
## cond.CategoryNon-Gap:PC1 0.6601346 2.2014151 1.2055000
## cond.CategoryOther Gap:PC1 1.2296463 2.5975110 1.7871821
## Location.cond.Std.Dev.(Intercept) 1.1056405 2.7882306 1.3783666
## zi.(Intercept) 0.1102974 1.5340883 0.4113465
## zi.CategoryNon-Gap 1.4069059 23.5968903 5.7618230
## zi.CategoryOther Gap 0.2640439 3.7872731 1.0000033
## Location.zi.Std.Dev.(Intercept) 1.5832860 14.1600376 3.0149229
I think I’ve figured out how to correctly interpret the odds ratios and intercept scores
#Intercept
#Non-Gaps
33.715419*0.1132062 %>% round(2)#lower
## [1] 3.708696
33.715419*0.2325890 %>% round(2) #estimate
## [1] 7.754546
33.715419*0.4778683 %>% round(2)#upper
## [1] 16.1834
#Other-Gaps
33.7154190*0.3387368 %>% round(2)#lower
## [1] 11.46324
33.7154190*0.5355936 %>% round(2)#estimate
## [1] 18.20633
33.7154190*0.8468538 %>% round(2)#upper
## [1] 28.65811
#PC1
#Non-Gaps
0.7128030*0.6601346 %>% round(2)#lower
## [1] 0.47045
0.7128030*1.2055000 %>% round(2)#estimate
## [1] 0.8624916
0.7128030*2.2014151 %>% round(2)#upper
## [1] 1.568167
#Other-Gaps
0.7128030*1.2296463 %>% round(2)#lower
## [1] 0.8767477
0.7128030*1.7871821 %>% round(2)#estimate
## [1] 1.275917
0.7128030*2.5975110 %>% round(2)#upper
## [1] 1.853288
emtrends extracts the slopes which are expressed relative to EAB gaps - I then use exp to “unlink” the predictions.
These values are really just PC1 * (Category:PC1)
I think the levels in here get mixed up somehow because the other-gap and non-gap categories come out mixed up - I’m going to use the values from above but I know where they’re coming from onw.
emtrends(shrubs.hurdle5, pairwise ~ Category, var = "PC1")
## $emtrends
## Category PC1.trend SE df lower.CL upper.CL
## EAB Gap -0.339 0.158 60 -0.654 -0.0228
## Other Gap -0.152 0.310 60 -0.771 0.4681
## Non-Gap 0.242 0.177 60 -0.111 0.5953
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## EAB Gap - Other Gap -0.187 0.307 60 -0.608 0.8162
## EAB Gap - Non-Gap -0.581 0.191 60 -3.044 0.0096
## Other Gap - Non-Gap -0.394 0.303 60 -1.298 0.4017
##
## P value adjustment: tukey method for comparing a family of 3 estimates
exp(-0.339) %>% round(3)
## [1] 0.712
exp(-0.152) %>% round(3)
## [1] 0.859
exp(0.242) %>% round(3)
## [1] 1.274
str(shrubs.EI)
## 'data.frame': 72 obs. of 20 variables:
## $ Location : chr "Code Farm" "Code Farm" "Code Farm" "Code Farm" ...
## $ Category : Factor w/ 3 levels "EAB Gap","Other Gap",..: 1 1 1 1 3 3 3 3 2 2 ...
## $ ID : int 1 2 3 4 1 2 3 4 1 2 ...
## $ richness : int 2 3 2 3 3 2 3 2 3 1 ...
## $ H : num 0.5 0.965 0.693 1.099 0.736 ...
## $ percentNN : num 0 0 0 0 12.5 ...
## $ percentB : num 0 0 0 0 12.5 ...
## $ Buckthorn : int 0 0 0 0 10 0 0 0 10 0 ...
## $ Edge.Distance : num 5.59 61.89 59.79 14.57 85.46 ...
## $ Fragment.Size : int 37 37 37 37 37 37 37 37 37 37 ...
## $ Percent.Tree : int 50 90 80 30 70 50 110 100 60 40 ...
## $ Management : chr "Private" "Private" "Private" "Private" ...
## $ Buckthorn_Management: int 0 0 0 0 0 0 0 0 0 0 ...
## $ Shade.Tolerance.Site: num 4.42 4.42 4.42 4.42 4.42 ...
## $ Gap.Area.Site : num 136 136 136 136 136 ...
## $ PC1 : num -1.383 -0.828 -0.826 -0.899 -1.083 ...
## $ PC2 : num 0.1294 -0.0277 0.2463 0.8533 -0.1041 ...
## $ PC1.s : num [1:72, 1] -1.026 -0.614 -0.613 -0.666 -0.803 ...
## ..- attr(*, "scaled:center")= num -1.89e-16
## ..- attr(*, "scaled:scale")= num 1.35
## $ res : num -0.582 -0.581 -0.581 -0.581 0.272 ...
## $ fit : num 71.26 47.5 47.44 50.03 6.76 ...
shrubs.summary <- shrubs.EI %>% select(Location, Category, ID, Buckthorn, Edge.Distance, Fragment.Size, Percent.Tree, Shade.Tolerance.Site, PC1) %>%
group_by(Location, Category) %>%
summarize(Buckthorn = mean(Buckthorn), Edge.Distance = mean(Edge.Distance), Fragment.Size = mean(Fragment.Size),
Percent.Tree = mean(Percent.Tree), Shade.Tolerance.Site = mean(Shade.Tolerance.Site), PC1 = mean(PC1)) %>%
mutate_if(is.numeric, round, digits = 1)
## `mutate_if()` ignored the following grouping variables:
## Column `Location`
shrubs.summary
## # A tibble: 18 x 8
## # Groups: Location [6]
## Location Category Buckthorn Edge.Distance Fragment.Size Percent.Tree
## <chr> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Code Fa… EAB Gap 0 35.5 37 62.5
## 2 Code Fa… Other G… 5 102. 37 60
## 3 Code Fa… Non-Gap 2.5 75.9 37 82.5
## 4 Field R… EAB Gap 55 62 11 92.5
## 5 Field R… Other G… 7.5 60.3 11 42.5
## 6 Field R… Non-Gap 5 57.4 11 52.5
## 7 Five Po… EAB Gap 30 192. 174 82.5
## 8 Five Po… Other G… 7.5 187. 174 60
## 9 Five Po… Non-Gap 2.5 201. 174 60
## 10 Meadowl… EAB Gap 7.5 257. 324 62.5
## 11 Meadowl… Other G… 2.5 256. 324 85
## 12 Meadowl… Non-Gap 0 248 324 67.5
## 13 Medway … EAB Gap 17.5 110. 230 100
## 14 Medway … Other G… 27.5 195. 230 70
## 15 Medway … Non-Gap 0 136. 230 70
## 16 Westmin… EAB Gap 47.5 55.8 199 115
## 17 Westmin… Other G… 37.5 90.3 199 115
## 18 Westmin… Non-Gap 10 80.7 199 105
## # … with 2 more variables: Shade.Tolerance.Site <dbl>, PC1 <dbl>
shrubs.summary$Location <- revalue(shrubs.summary$Location, c("Field Research Station"="Private 1"))
shrubs.summary$Location <- revalue(shrubs.summary$Location, c("Code Farm"="Private 2"))
shrubs.summary$Location <- revalue(shrubs.summary$Location, c("Westminster Ponds"="Public 1"))
shrubs.summary$Location <- revalue(shrubs.summary$Location, c("Five Points Forest"="Private 3"))
shrubs.summary$Location <- revalue(shrubs.summary$Location, c("Medway Valley"="Public 2"))
shrubs.summary$Location <- revalue(shrubs.summary$Location, c("Meadowlily Woods"="Public 3"))
#write.csv(shrubs.summary, file = "shrubs_summary.csv")
gap.area.summary <- EI5 %>% #summarize gap area
filter(Category != "Non-Gap") %>%
group_by(Location, Category) %>%
summarize(Gap.Area.Site = mean(Gap.Area) %>% round(1))
gap.area.summary$Location <- revalue(gap.area.summary$Location, c("Field Research Station"="Private 1"))
gap.area.summary$Location <- revalue(gap.area.summary$Location, c("Code Farm"="Private 2"))
gap.area.summary$Location <- revalue(gap.area.summary$Location, c("Westminster Ponds"="Public 1"))
gap.area.summary$Location <- revalue(gap.area.summary$Location, c("Five Points Forest"="Private 3"))
gap.area.summary$Location <- revalue(gap.area.summary$Location, c("Medway Valley"="Public 2"))
gap.area.summary$Location <- revalue(gap.area.summary$Location, c("Meadowlily Woods"="Public 3"))
#write.csv(gap.area.summary, file = "gaps_summary.csv")
As per reviewer comments - influence of this station
Remove FRS from data
shrubs.EI.wo <- shrubs.EI %>% filter(Location != "Field Research Station") %>%
select(Location, Category, ID, Buckthorn, Management, Gap.Area.Site, PC1)
str(shrubs.EI.wo)
## 'data.frame': 60 obs. of 7 variables:
## $ Location : chr "Code Farm" "Code Farm" "Code Farm" "Code Farm" ...
## $ Category : Factor w/ 3 levels "EAB Gap","Other Gap",..: 1 1 1 1 3 3 3 3 2 2 ...
## $ ID : int 1 2 3 4 1 2 3 4 1 2 ...
## $ Buckthorn : int 0 0 0 0 10 0 0 0 10 0 ...
## $ Management : chr "Private" "Private" "Private" "Private" ...
## $ Gap.Area.Site: num 136 136 136 136 136 ...
## $ PC1 : num -1.383 -0.828 -0.826 -0.899 -1.083 ...
Run model without FRS data
shrubs.hurdle5.wo <- glmmTMB(Buckthorn ~ Category * PC1 + (1 | Location),
ziformula = ~ Category + (1 | Location), family = "truncated_nbinom2", data = shrubs.EI.wo)
exp(fixef(shrubs.hurdle5)$cond) # Intercept expected cover when present, others are multiplicative changes in cover conditional on presence
## (Intercept) CategoryNon-Gap CategoryOther Gap
## 33.7154190 0.2325890 0.5355936
## PC1 CategoryNon-Gap:PC1 CategoryOther Gap:PC1
## 0.7128030 1.2055000 1.7871821
exp(fixef(shrubs.hurdle5.wo)$cond)
## (Intercept) CategoryOther Gap CategoryNon-Gap
## 49.2315112 0.3845417 0.1975805
## PC1 CategoryOther Gap:PC1 CategoryNon-Gap:PC1
## 0.4407283 3.7069936 2.2689714
Remove FRS:
Therefore, greater EAB & EI effect when we remove FRS
#Non-gaps
33.7154190 * 0.2325890
## [1] 7.841836
49.2315112 * 0.1975805
## [1] 9.727187
#Other gaps
33.7154190 * 0.5355936
## [1] 18.05776
49.2315112 * 0.3845417
## [1] 18.93157
exp(fixef(shrubs.hurdle5)$zi) # Intercept is odds of zero, others are odds ratios
## (Intercept) CategoryNon-Gap CategoryOther Gap
## 0.4113465 5.7618230 1.0000033
exp(fixef(shrubs.hurdle5.wo)$zi) # Intercept is odds of zero, others are odds ratios
## (Intercept) CategoryOther Gap CategoryNon-Gap
## 0.5730672 0.7672057 4.9923899
Remove FRS:
exp(confint(shrubs.hurdle5))
## 2.5 % 97.5 % Estimate
## cond.(Intercept) 22.3264815 50.9139553 33.7154190
## cond.CategoryNon-Gap 0.1132062 0.4778683 0.2325890
## cond.CategoryOther Gap 0.3387368 0.8468538 0.5355936
## cond.PC1 0.5231179 0.9712689 0.7128030
## cond.CategoryNon-Gap:PC1 0.6601346 2.2014151 1.2055000
## cond.CategoryOther Gap:PC1 1.2296463 2.5975110 1.7871821
## Location.cond.Std.Dev.(Intercept) 1.1056405 2.7882306 1.3783666
## zi.(Intercept) 0.1102974 1.5340883 0.4113465
## zi.CategoryNon-Gap 1.4069059 23.5968903 5.7618230
## zi.CategoryOther Gap 0.2640439 3.7872731 1.0000033
## Location.zi.Std.Dev.(Intercept) 1.5832860 14.1600376 3.0149229
exp(confint(shrubs.hurdle5.wo))
## 2.5 % 97.5 % Estimate
## cond.(Intercept) 28.52113323 84.9805538 49.2315112
## cond.CategoryOther Gap 0.18017215 0.8207279 0.3845417
## cond.CategoryNon-Gap 0.08133583 0.4799613 0.1975805
## cond.PC1 0.24618336 0.7890112 0.4407283
## cond.CategoryOther Gap:PC1 1.45022580 9.4756291 3.7069936
## cond.CategoryNon-Gap:PC1 0.53113419 9.6929014 2.2689714
## Location.cond.Std.Dev.(Intercept) 1.00000000 Inf 1.0000605
## zi.(Intercept) 0.13614110 2.4122474 0.5730672
## zi.CategoryOther Gap 0.18368121 3.2044900 0.7672057
## zi.CategoryNon-Gap 1.06280095 23.4511994 4.9923899
## Location.zi.Std.Dev.(Intercept) 1.55107946 19.6450909 3.1370843
Remove FRS:
pred.wo <- data.frame(PC1 = c(seq(EAB.min, EAB.max, length.out = 100),
seq(other.min, other.max, length.out = 100),
seq(non.min, non.max, length.out = 100)),
Category = c(rep("EAB Gap", 100), rep("Other Gap", 100), rep("Non-Gap", 100)),
Gap.Area.Site = rep(mean(shrubs.EI$Gap.Area.Site),300), # Added gap area
Management = c(rep("Private", 150), rep("Public", 150)),
Location = rep(NA, 300)) # This NA gives general predictions instead of location specific predictions
pred.wo$fit <- predict(shrubs.hurdle5.wo, newdata = pred, se.fit = TRUE, type = "link")$fit
pred.wo$se <- predict(shrubs.hurdle5.wo, newdata = pred, se.fit = TRUE, type = "link")$se.fit
str(pred.wo)
## 'data.frame': 300 obs. of 7 variables:
## $ PC1 : num -2.28 -2.23 -2.19 -2.14 -2.1 ...
## $ Category : chr "EAB Gap" "EAB Gap" "EAB Gap" "EAB Gap" ...
## $ Gap.Area.Site: num 139 139 139 139 139 ...
## $ Management : chr "Private" "Private" "Private" "Private" ...
## $ Location : logi NA NA NA NA NA NA ...
## $ fit : num 5.76 5.73 5.69 5.65 5.62 ...
## $ se : num 0.886 0.873 0.86 0.847 0.834 ...
pred.upr <- pred.wo$fit + 1.96 * pred$se
pred.lwr <- pred.wo$fit - 1.96 * pred$se
pred.wo$upr <- ifelse(exp(pred.upr) < 100, exp(pred.upr), 100)
pred.wo$lwr <- ifelse(exp(pred.lwr) > 0, exp(pred.lwr), 0)
pred.wo$mean <- exp(pred.wo$fit)
shrubs.EI.wo$Category = factor(shrubs.EI.wo$Category, levels=c("EAB Gap", "Other Gap", "Non-Gap"))
pred.wo$Category = factor(pred.wo$Category, levels=c("EAB Gap", "Other Gap", "Non-Gap"))
fig.EI.wo <- ggplot(pred.wo, aes(x = PC1, y = mean)) +
geom_ribbon(aes(x = PC1, ymin = lwr, ymax = upr, fill = Category), alpha = 0.15) +
geom_point(data = shrubs.EI.wo, aes(x = PC1, y = Buckthorn, col = Category), alpha = 0.6) +
geom_line(aes(col = Category), size = 1, alpha = 0.8) +
labs(x = "Ecological Integrity",
y = "Buckthorn (%)",
col = "Category") +
scale_y_continuous(limits = c(0, 100)) +
scale_colour_manual(values = mypalette2) +
scale_fill_manual(values = mypalette2) +
theme_classic() +
theme(legend.position = "none")
fig.EI
fig.EI.wo
## Warning: Removed 32 row(s) containing missing values (geom_path).
#ggsave('figures/fig.EI.noFRS.jpeg',fig.EI.wo , units = 'cm', width = 14, height = 10)
Sys.time()
## [1] "2020-05-22 15:05:38 PDT"
git2r::repository()
## Local: master /Users/JenBaron/Documents/UWO/UWO 4th Year/Publication/Analysis/EAB_Manuscript
## Remote: master @ origin (https://github.com/JenBaron/EAB_Manuscript.git)
## Head: [3b334fb] 2020-05-22: figure with legend
sessionInfo()
## R version 4.0.0 (2020-04-24)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] emmeans_1.4.6 DHARMa_0.3.0 gridExtra_2.3 ggrepel_0.8.2
## [5] RColorBrewer_1.1-2 glmmTMB_1.0.1 lme4_1.1-23 Matrix_1.2-18
## [9] nlme_3.1-147 forcats_0.5.0 stringr_1.4.0 dplyr_0.8.5
## [13] purrr_0.3.4 readr_1.3.1 tidyr_1.0.2 tibble_3.0.1
## [17] ggplot2_3.3.0 tidyverse_1.3.0 plyr_1.8.6
##
## loaded via a namespace (and not attached):
## [1] fs_1.4.1 lubridate_1.7.8 doParallel_1.0.15
## [4] httr_1.4.1 numDeriv_2016.8-1.1 tools_4.0.0
## [7] TMB_1.7.16 backports_1.1.6 utf8_1.1.4
## [10] R6_2.4.1 mgcv_1.8-31 DBI_1.1.0
## [13] colorspace_1.4-1 withr_2.2.0 tidyselect_1.0.0
## [16] git2r_0.26.1 compiler_4.0.0 cli_2.0.2
## [19] rvest_0.3.5 ROI_0.3-3 xml2_1.3.2
## [22] spaMM_3.2.0 labeling_0.3 slam_0.1-47
## [25] scales_1.1.0 mvtnorm_1.1-0 pbapply_1.4-2
## [28] proxy_0.4-24 digest_0.6.25 minqa_1.2.4
## [31] rmarkdown_2.1 pkgconfig_2.0.3 htmltools_0.4.0
## [34] fastmap_1.0.1 dbplyr_1.4.3 rlang_0.4.5
## [37] readxl_1.3.1 rstudioapi_0.11 shiny_1.4.0.2
## [40] farver_2.0.3 generics_0.0.2 jsonlite_1.6.1
## [43] magrittr_1.5 Rcpp_1.0.4.6 munsell_0.5.0
## [46] fansi_0.4.1 lifecycle_0.2.0 Rglpk_0.6-4
## [49] stringi_1.4.6 yaml_2.2.1 MASS_7.3-51.6
## [52] grid_4.0.0 promises_1.1.0 qgam_1.3.2
## [55] parallel_4.0.0 crayon_1.3.4 lattice_0.20-41
## [58] haven_2.2.0 splines_4.0.0 hms_0.5.3
## [61] knitr_1.28 pillar_1.4.3 boot_1.3-25
## [64] estimability_1.3 codetools_0.2-16 reprex_0.3.0
## [67] glue_1.4.0 evaluate_0.14 modelr_0.1.6
## [70] httpuv_1.5.2 vctrs_0.2.4 nloptr_1.2.2.1
## [73] foreach_1.5.0 cellranger_1.1.0 gtable_0.3.0
## [76] assertthat_0.2.1 xfun_0.13 mime_0.9
## [79] xtable_1.8-4 broom_0.5.6 later_1.0.0
## [82] coda_0.19-3 iterators_1.0.12 registry_0.5-1
## [85] ROI.plugin.glpk_0.3-0 statmod_1.4.34 ellipsis_0.3.0
## [88] gap_1.2.2